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System for Estimating Azimuthal Variations in Seismic Data 
BACKGROUND OF THE INVENTION 

1. Field of the Invention 

This invention is related to seismic data processing. More specifically, the 
invention is related to a system for processing seismic data to detect azimuthal velocity 
variations. 

2. Descriptionof Related Art 

Seismic surveys are routinely used in the search for oil and gas reservoirs in the 
earth's subsurface. Seismic surveys are performed by imparting acoustic energy into the 
earth, either at a land surface or in a marine environment, and then detecting the reflected 
and refi*acted acoustic energy. The delay time between the imparting of the acoustic energy 
wave at the source location and detection of the same wave at a receiver location indicates 
the depth of reflecting geological interfaces. 

Until recently, only two-dimensional ("2D") seismic surveys were conducted, with 
the seismic source locations being coUinear with a line of receivers. Recent advances in 
technology have enabled three-dimensional ("3D") seismic survey data to be gathered and 
analyzed. Typically, in 3D surveys arrays of seismic receivers are deployed which receive 
reflected acoustic energy imparted at varying locations that may be specifically selected to 
provide a rich assortment of azimuths for common midpoints. 

A technique firequently used in seismic survey analysis is AVO analysis, which is 
the amplitude variation with offset, and is also referred to herein as amplitude variation 
with incidence angle. According to the AVO approach, attributes of a subsxirface interface 
are determined both firom the normal-incidence amplitude of reflected seismic energy, and 
also fi-om the dependence of the detected seismic reflections on the angle of incidence of 
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the seismic energy at a subsurface reflecting interface relative to the vertical. In 
conventional AVO analysis, multiple seismic traces having a common reflection point, 
commonly referred to as a common mid point or common depth point(CMP or CDP) ' 
gather, are collected. From the CMP (or CDP) gather, one may derive the amplitude R of a 
reflected seismic wave from an interface (i.e., the "target horizon") as a function of the 
angle of incidence 6 from the normal according to the following relationship: 

In this case, the coefficient A is the zero-offset response (also referred to as the AVO 
intercept), while the coefficient B is referred to as the AVO slope, or gradient, as it is 
representative of the rate of change of amplitude with the square of the angle of incidence 
Analysis of the AVO slope and intercept can provide indicators of interesting formations, 
from an oil and gas exploration standpoint. For example, variations in the A and B values 
from a theoretical A-versus-B trend line for the expected stratigraphic sequences can 
indicate the location of hydrocarbon reserves. 

While simple models of subsurface geology assume azimuthal isotropy m the 
propagation of acoustic energy it has been observed that azimuthal anisotropy is in fact 
present in many survey regions, such that the velocity of acoustic energy depends upon the 
azimuth of the source-receiver path. If azimuthal anisotropy is present, the conventional 
normal moveout correction may not adequately align the seismic traces in the gather, which 
can result in degraded AVO analysis. 

Normal moveout correction of the seismic data, both for offset-dependent delays 
and also for azimuthal anisotropy caused by the overburden, is therefore typically 
performed in producing stacked traces of improved signal-to-noise ratio for use in a 3D 
seismic survey. For example. U. S. Patent No. 5,532,978 describes a method of deriving 
and applying azimuthal anisotropy corrections to seismic survey signals. 
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The detection of a preferred azimuthal direction at a reflecting interface can also 
provide important information regarding geological features. For example, a preferred 
azimuthal reflection direction can indicate the presence of aligned vertical fractures. For 
moderately far offsets (25°- 35° incidence angles), the P-wave traveling in the plane wave 
5 parallel to aligned vertical fractures has a higher velocity than the P-wave traveling in the 
plane perpendicular to the fractures. 

Traditionally, azimuthal velocity analysis has been performed using azimuth- 
sectored supergathers and picking semblance maxima at various azimuths. This reduces 
the problem to a series of 2-D solutions, rather than solving the complete 3-D solution. In 
10 some cases as few as two sectors may be chosen, perpendicular and parallel to the 

(average) principal axes of the azimuthal anisotropy. If more than two sectors are used, an 
ellipse is fitted to the picked velocities to give fast and slow velocity magnitudes and the 
azimuth of the fast velocity. These procedures suffer from several drawbacks: 

Picking semblance, by hand, on azimuth sectored data is processor/interpreter 
1 5 dependent and extremely time consuming. 

Semblance works well for data which do not show amplitude variation with offset 
("AVO"), however, if the data contain significant AVO, particularly if there is a polarity 
reversal, semblance can fail In this case automatic picking of semblance maxima will be 
erroneous. 

20 If the subsurface has azimuthal velocity variation ("AW") then this will appear as 

an offset-dependent static viewed on offset-sorted CMP gathers. This will reduce the 
effectiveness of any surface consistent statics solution, thus the azimuth-sectored 
supergathers will most likely be contaminated with statics. This will significantly degrade 
the semblance analysis and may result in several semblance maxima. 
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The semblance is based on giving the greatest stack power. However, for AW 
analysis it is the actual subsurface velocity that is of interest, not simply the velocity that 
gives the best stack. For instance, if a higher amplitude occurs at a particular azimuth 
within the sector, then the velocity at that azimuth will be picked. In addition, if those high 
5 amplitudes are at the mid to near offsets and are contaminated with residual statics then a 
completely erroneous velocity could give the highest semblance. 

Sectoring and partial stacking of the data means that it is extremely difficult to 
obtain error estimates. Not only is it difficult to attribute a picking error from picking 
semblance, but errors due to the acquisition geometry are not represented. In any analysis 
of this type it is important to compute the errors associated with the obtained results. For 
instance a weighted least squares approach has been used to compute the errors in a 
technique for inverting azimuthal variation of amplitude for shear wave data. It has also 
been observed that the reliability of the amplitude variation v^th azimuth analysis has been 
assessed by looking for an absence of the acquisition geometry being mirrored in the 
anisotropy maps. 

It should be noted that the description of the invention which follows should not be 
construed as limiting the invention to the examples and preferred embodiments shown and 
described. Those skilled in the art to which this invention pertains will be able to devise 
variations of this invention within the scope of the appended claims. 

20 SUMMARY OF THE INVENTION 

The invention comprises a system for processing seismic data to estimate time shift 
resulting from velocity anisotropy in the earth's subsurface. A gather of seismic data traces 
is formed and selected seismic data traces included in said gather within selected time 
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windows are cross-correlated to estimate the time shift in the seismic data traces included 
in said gather resulting from velocity anisotropy in the earth's subsurface. 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 shows a simplified portion of a 3D seismic source-receiver layout for a 3D 
5 seismic survey. 

FIG. 2A shows representative seismic data traces from a common midpoint gather 
prior to application of normal moveout adjustment. 

FIG. 2B shows the seismic data traces shown in FIG. after application of normal 
moveout adjustment. 
10 FIG. 3 is a flow chart showing an embodiment of the invention. 

FIG. 4 is a flow chart showing a further embodiment of the invention. 
FIG. 5 is a flow chart showing a further embodiment of the invention. 
FIG. 6 is a flow chart showing a further embodiment of the invention.. 
FIG. 7 is a flow chart showing a further embodiment of the invention. 
15 FIG. 8 is a diagram illustrating a spatial relationship of the invention. 

FIG. 9 shows a computer system for carrying out the invention. 

DESCRIPTION OF PREFERRED EMBODIMENT 

FIG. 1 shows a simplified portion of a 3D seismic source-receiver layout for a 3D 
20 seismic survey. Shown in FIG. 1 is a portion of a receiver array 12, positioned on the 

earth's surface, comprising three columns of receivers, with each column including eight 
receivers. Source array 14 includes a group of source locations. Typically, a seismic 
source is moved along the surface of the earth, and the source is activated at specific source 
locations in a sequence. The acoustic energy imparted at each source location travels 
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U^ugh *e canh and. after reflection torn subsurface geologic interfaces, is detected by 
each receiver in the receiver array. 

Also shown in FIG. 1 is an example ofa common midpoint MP. The common 
midpointshowniscommontoseveralsourcereceiverpathsinthissurvey. FIG. 1 
illustrates MP as being a midpoint between source locadon SU and recerver loca«on R38. 
^.urce location S21 and receiver location R28; source iocaUon S38 and receiver locatron 
RU. and source locaUon S36 and receiver location R13. It will be apparent to those of 
ordmary skiU in .he art that location MP could be the common midpoint for a large number 
of sotmx-receiver pairs, and also that there are a large number of common midpomts 
betweenothersourcer^eiverp^. Each ofthese source-receiver paths may be of a 
diffe^nt length(offsc.)and different direction(azin,u.h).Alth„ughcommonm.d^^^ 

^ are typically ^ in implementing the present invemion. if «.e 

dipping rather than flat, additional processingknown to .hoseofordi™ryslt.U.nthe art 

i;jperfor»edonfl,eseismicdatatodeve.opc„mntonreflectionpointgathers.owh.ch 
the present invention may be applied. 

AS is well known in the art, normal moveont corrections are typieaUy made to 
.races inacommon midpoint gather to correct for the additional delay .im. f«lo.^« 
ome.traces.sothatthetravcl.ime„feachseismicsignaliseffective,ynorn« 

^ro^flse. travel time. In sin^tions where the earth exhibits no azmtutha. amsotropy, the 
^uthalvariationwiUnotintroduceavariaUoninthedataset. However, when 
a^uthal aniso^opy is present in seismic survey signals, the standard normal moveout 
correction may no. adequately correct for variations in delay times for traces from sour^ 
receiver pairs of different directions. 

variations in the seismic datatraces will appear to be variations in ampliti.de as a 
«^onof the — variation, wheninrcalitytitevariationsinti^sei^tic^^^^ 
3rearesuI.ofaximuti«llydependenttimeshifls.Avelocityvariationw,ti.azm,uti.of„nly 



a few percent can cause ten milliseconds or more of time difference in the location of an 
event in a seismic data trace. Accordingly, velocity variations below the resolution of 
conventional semblance based velocity analysis can distort the amplitude variation with 
incidence angle (AVO) and the amplitude variation with azimuth ("AVOA"), resulting in 
5 an incorrect computation. 

In addition to the standard normal moveout correction, typically, for 3-D land 
acquisition, deconvolution, refraction statics, and two passes of surface consistent residual 
statics processes are applied to the data. As a preliminary to the performance of the 
10 invention as disclosed herein, noise reduction, trace scaling and any other single-trace 
process may be applied, however, the application of multi-trace processing should be 
avoided. 

FLATTEN EVENTS WITHIN A MOVING TIME WINDOW 

If the velocity of the seismic signals in the subsurface vary with the azimuthal 

15 direction in which the seismic signals travel, the conventional normal moveout correction 
will not align the traces properly, and in accordance with a preferred embodiment of the 
invention, the following process may be utilized to achieve proper trace alignment. FIG. 
2A and FIG. 2B show representative traces from a common midpoint gather. An initial 
velocity profile for the seismic traces may be generated in a conventional manner and 

20 normal moveout applied to the traces. FIG. 2A shows the traces prior to applying normal 
moveout adjustment to the traces and FIG. 2B shows these same traces after application of 
the normal moveout adjustment. As outlined in FIG. 3, and with reference to FIG. 2B, in a 
preferred embodiment of the invention a first time window is selected in step 50. Time 
windows utilized for performing the invention may typically be in the range of 100 to 300 

25 milliseconds, however, the time window utilized may vary in accordance v^th the 

judgment of the processor. A first selected time window for the data shown in FIG. 2B 
may extend from 500 milliseconds to 700 milliseconds. In step 52, a spatial window 
comprising a number of spatially related traces is selected and the traces in this spatial 
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window are summed together within the selected time window to create a "pilot" trace. In 
step 54, one trace is selected from the pilot trace, which will be referred to as the "input" 
trace. In step 56 the summed pilot trace is cross-correlated with the input trace within the 
selected time window and in step 58 the time shift between the input trace and pilot trace 
5 which yields the maximum correlation is determined and recorded. 

The spatial window is then moved across the CMP gather within the selected time 
window by one trace and a new pilot trace is generated for this new spatial window. The 
new pilot trace is then correlated with the new input trace within the time window. This 
process is continued for the selected time window, with each successive data trace within 

10 the gather being designated as the "input trace" and cross-correlated with a pilot trace 

which comprises a plurality of nearby traces, to complete a time window trace correlation 
sequence. In this way the pilot trace represents the local phase and amplitude 
characteristics of the data for each "input" trace. Accordingly, a decision is made in step 
60 as to whether a time window trace correlation sequence has been completed. If the 

15 answer is no, steps 52, 54, 56, 58 and 60 are repeated for a successive input trace. If the 
answer is yes, then in step 62 a decision is made as to whether the cross-correlation 
sequence just completed is the first performed cross-correlation process for the time 
window. In one implementation of the invention, if the answer is yes, the calculated time 
shift that provided maximum cross-correlation between the input trace and the pilot trace 

20 for each trace is applied to each trace in step 64 and steps 52, 54, 56, 58, 60 and 62 are 
repeated. If the answer in step 62 is no, then in step 66 a decision is made as to whether 
any reflection event in the time window is substantially aligned across the traces in the 
gather. If the answer is no steps 64, 52, 54, 56, 58, 60, 62 and 66 are repeated. This 
decision in step 66 is normally based on whether or not the additional time shifts calculated 

25 for the input traces in the just completed time window trace correlation sequence are 

significant. If the answer in step 66 is yes, then in step 68, the total computed time shifts 
for each trace are stored and the traces are returned to their form they were in prior to 
beginning the process outlined in FIG. 3. Typically, only two iterations of the process 
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described in steps 52, 54, 56, 58, 60, 62 and 64 are performed, but further iterations may be 
performed if, in the judgment of the processor, data quality may be improved by further 
iterations. 

A decision is made in step 70 as to whether ail time windows of interest in the 
5 seismic data gather have been selected, and if the answer in no a new time window is 
selected in step 50 and the cross-correlation procedure described above with respect to 
steps 52, 54, 56, 58, 60, 62 64, 66, 68 and 70 is repeated for all time v^ndows of interest. 
The successive time windows selected may occupy successive time positions on the 
seismic data traces or the time windows may overlap, depending on the quality of the data. 

10 Following completion of the cross-correlation procedure for all time windows and 

all traces within each time window, in step 72, the amount of time shift which achieved 
the maximum correlation for each trace within each time window is applied to each trace at 
the center point within each time window, and, in step 74, time shifts for the remainder of 
the data traces are interpolated between these center points. 

1 5 Each pilot trace may comprise, for example, eleven traces. The trace in the center 

of the selected spatial window, i.e. the sixth trace, may be designated as the "input" trace 
and cross-correlated with the pilot trace to obtain a time shift. The spatial trace v^ndow is 
then moved one trace across the gather and a new pilot trace formed. Again, the trace in 
the center of this window, i.e. the next trace in the gather, is designated as the "input" trace 

20 and cross-correlated with the pilot trace to obtain a time shift for that trace. At the edges of 
the gather, in this example the first through the sixth traces, the spatial window comprising 
the pilot trace may be shortened so that, if the first trace is the "input" trace, the first 
through sixth traces are stacked to form the pilot trace, and for the second trace, the first 
through the seventh traces are stacked and so on, until the full number of traces desired in 

25 the spatial window is reached (in this example, eleven). The number of traces selected to 
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form the pilot trace may be selected on the basis of the magnitude of amplitude variation 
with offset in the data and the magnitude of noise such as muhiple contamination. 

In performing the cross-correlation procedure, amplitude and phase variations 
with offset, including the case where an event reverses polarity at some offset are taken 
into accoimt. The pilot trace represents the "local" data characteristics. If an event 
reverses polarity at far offsets then the pilot trace at near offsets should not include the far 
offset traces. Similarly, the pilot trace at far offsets will not include traces at near offsets. 

Other trace attributes, including absolute values, RMS values, or trace envelope 
may be used for performing the cross-correlation in addition to the raw trace reflection 
amplitude. 

At such time as the time shifts have been applied to the seismic data traces in the 
gather, AVO analysis as well as AVOA analysis, such as discussed herein with reference to 
FIG. 4, may be performed on the adjusted traces. 

COMPUTING THE AMPLITUDE VARIATION WITH INCIDENCE ANGLE (AVO), 
AND THE AMPLITUDE VARIATION WITH AZIMUTH (AVOA) 

It is known to those of ordinary skill in the art that amplitude variation with 
incidence angle (also referred to as amplitude variation with offset), as well as the 
amplitude variation with azimuth for a reflection from a horizontal transverse isotropic 
layer of the earth which is overlain with an isotropic overburden can be approximated by 
the following equation: 

R{dJ)= I+G,Sm^(0)+G2Sm\0)Cos\^-fi) (Eq. 1) 
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where: R(0,^) is the reflection coefficient as a function of 0 , the incidence angle of the 
seismic energy at a subsurface reflecting interface relative to the vertical, and ^ , 

the receiver azimuth with respect to a predefined zero azimuth direction (for 
example, true north); 

/ is the P-wave impedance contrast between the subsurface layers from which the 
signal is reflected; 

Gj is the isotropic AVO gradient; 

G2 is the azimuthal or anisotropic term; and 

(with reference to FIG. 8 ) is the angle between the predefined zero azimuth 
direction (such as true north) and the maximum AVO gradient direction. 

Gj and are given by: 



where: A/7 , A and A are the change in density, the change in P-wave velocity, and 
the change in S-wave velocity, respectively, 



G. = 




(Eq. 2) 



and 




(Eq. 3) 
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p , and are the average density, the average P-vy^ave velocity and the 
average S-wave velocity respectively, 

g = , the average P-wave velocity divided by the average S-wave velocity, 



A S^^^ is the change in S across the reflecting boundary, and 
A / is the change in the shear wave splitting parameter / across the reflecting 
boundary, where: 



r = 



c - c 

^1212 ^3232 
2C3232 



It is known to those of ordinary skill in the art that for a linearly elastic material, each 
component of stress is linearly dependent on every component of strain e^^i , where / , 

j , k and / are directional indices that may assume values of 1, 2 or 3. The stress-strain 
dependency is given by Hooke's Law: 

where C^.^ is the elastic modulus tensor and completely characterizes the elasticity of the 
medium. The relationship between S and the elastic modulus tensor is given by: 

^(v) _ (Q133 ~ Q232) " (3333 "Q232) 
2^333(^333 " Q232) 

However, without knowing >^ , Eq. 1, cannot be solved using a least squares 
approach. However, Eq. 1 can be rewritten as: 

R{0J)^ I+{g, + G^Cos\(I) ' P)]Sin\e) (Eq.4) 



-12- 



AXG-001 



which can be rewritten as: 

sothat G, = G;andG2 = G2'-G=^• 
Utilizingtheequality■. 



then 



(Eq. 7) 



10 



It is known to those of ordinary skill in the art that: 
G," Co/ (j> - A) + Gi'SOiH^ - ^) = 

w W andfr which can be related back to the 
which is linear in the unknowns ;n,,M',j anaw.j.wn 

unknowns G;.G; and ^, as follows: 



(Eqll) 



T^us. combining E<,uationslwifl.E<,uatio„s4-ll.E,uaUonlcanbewri..enas: 
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with 



(Eq. 13) 



(Eq. 14) 



and 



>^ = ATAN 



Ml -^22)' + 4^2') 



(Eq.l5) 




Values of the reflection coefficient R(0,^) for specific values of the incidence 
angle 0 and the source-receiver azimuth ^ can be obtained firom the recorded seismic 

data for each reflection event by extracting the amplitudes of the seismic data traces as a 
function of offset and azimuth. With reference to FIG. 4, in step 80, values of the 
reflection coefficient R{0,^) and the source-receiver azimuth ^ are obtained fsom the 

seismic data that is being processed. To obtain the value of an incidence angle 6 , a 
smoothed version of the interval velocity is calculated in step 82 in a manner well known 
to those of ordinary skill in the art, and the RMS velocity is computed in step 84 fi-om the 
smoothed version of the interval velocity. In step 86, the values of the incidence angle, 0 , 
may then be determined, utilizing the following equation: 



0 = ASINi 




(Eq. 16) 
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where: 

X is the source to receiver offset; 

TJj is the zero offset two way travel time; 

V^^^ is the RMS velocity; and 

5 Fj^j is the interval velocity at the time of interest. 

In step 88, a least squares method is used to compute reflection coefficient as a function of 
azimuthal angle and incidence angle for the seismic traces comprising the CMP gather. 
Eq, 12 is solved in a straightforward least squares manner, known to those of ordinary skill 
in the art, for the unknowns, I .W^^.W^^.zsAW^^- Values for G, (the isotropic AVO 
10 gradient), (the azimuthal or anisotropic term) and /? may then be computed from 
values computed for , 1^22 ' ^12 • Accordingly, it is demonstrated that Eq. 1 is 
linear in / , Gj , G2 and the direction ^ . 



Note that as indicated above in Eqs. 2 and 3, the derived gradients G, and G2 are 

related to physical rock properties — , — , =, , A<^^^ and A/ . 

Vs P 



1 5 COMPUTING THE AZIMUTHAL VELOCITY VARIATION (AW) 

Because the process outlined in FIG. 3 determines the time shift in seismic traces 
associated with the azimuthal velocity variation, this determined time shift information 
may be used to compute the actual azimuthal velocity variation. Steps for computing the 
azimuthal variation of velocity are outlined in FIG. 5. In step 90 the total travel time T for 

20 each trace is computed by adding the time shifts determined in the process described with 
reference to FIG. 3 which achieved maximum correlation for each trace and the time shift 



-15- 



AXG-001 



obtained as a result of standard NMO correction to 7^ , the zero offset travel time. The 
following equation may then be utilized to solve for the azimuthal velocity variation: 



^ nmo \r J 



where: 

T = total travel time 

7^ = two way zero-ofifset traveltime 

X = offset 

^nmo (^^) = the azimuthally varying velocity as a function of the azimuth ^ , 

and 

1 



= W,,Cos\(l>) + 2W,^Cos{(l))Sin{(l>) + W^^Sin'iil)) (Eq. 18) 



Accordingly, the total traveltime may be written as: 

= + [w,,Cos^{(l>) + 2W,^Cos{(l>)Sin{(l>) + W^^Sin" {(l^)]x'' . (Eq. 19) 



In step 92, Eq. 19 may be solved by using a linear least squares method known to 
those of ordinary skill in the art, using the time shifts picked from the cross-correlation 
process described with reference to FIG. 3 which achieved maximum correlation. Eqs. 9, 

10 and 1 1 may then be used to obtain Gj* , Gj* and ^ . The fastest velocity and the 
slowest velocity are calculated from the calculated values of Gj* and Gj*. The fastest 



velocity is given by Vj-^^ - 
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1 

the slowest velocity is given by y^iow = i — r » and the azimuth of the slowest velocity is 

V^2 



given by J3 . 



Because the travel times are being fitted by the least squares solving of Eq. 19, the 
azimuth J] that is computed using Eq. 19 is the azimuth of the greater travel time. 

Accordingly, if the travel time is greater, the velocity is slower. The azimuthal velocity 
variation may then be computed in step 94 from the following relationship: 

1 1 ^ . . 1 



LT 2Cos\(l>-P)^—^Sin\(l>-P). (Eq.20) 

*^rmo \r) ^slow ^ fast 

The amount of time shift resulting from azimuthal anisotropy may be determined 
for each reflection event as a fiinction of azimuthal angle and the appropriate time shift 
10 may be applied to each trace to adjust for the azimuthal time shift. At such time as the time 
shifts have been applied to the seismic data traces in the gather, AVO analysis as well as 
AVOA analysis, such as discussed herein with reference to FIG. 4, may be performed on 
the adjusted traces. 

COMPUTING OF ERRORS ASSOCIATED WITH CALCULATION OF TIME SfflFT 
15 VARIATION, VELOCITY VARIATION AND AMPLITUDE VARL\TION WITH 
AZIMUTH 

Errors associated with the calculation of the time shift variation with azimuth, the 
velocity variation with azimuth, and the amplitude variation with azimuth may be 
estimated utilizing a least squares approach. In step 100, the least squares approach to 
20 estimating the errors associated with the calculation of the time shift variation with azimuth 
is formulated in matrix notation, and may be written: 

Ax^b (Eq.21) 
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Where b is a 1 x N matrix (i.e. a vector) containing tlie data (e.g. travel times or 

amplitudes), A is an M x N matrix of coefficients (e.g. Sin^ {9) ) and x is a 1 x M 

matrix (i.e. a vector) of the parameters to be solved for. For instance for Eq. 19: 
^1 A'l^cos^^i A'l^cos^i sin^, XySim^^^^ 



A = 



1 A'2COS^^2 -Sfj cos^2 sin^Jj ^zsin^^j 
1 Xlcos^<l>i A'j^ cos ^3 sin ^Jj Xlsin^^^ 



1 Xlcos^(l>„ 2Jr^os«J„sin^„ ^„'sin^<*J 



^0 



11 



and 



f f2 \ 



b = 
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Where 71 are the observed travel times for data with source-receiver offsets and 
azimuth . The matrix equation is equivalent to N simultaneous equations for M 

xmknowns. In the example shown, M is equal to four. For a least squares formulation, N, 
the number of data points, must be greater than M, the number of unknowns. There are 
5 various standard numerical methods, known to those of ordinary skill in the art, for 
computing the M unknowns. 

In step 102, the standard error in the unknowns is computed by taking the square 
root of the diagonals of the matrix E: 

E^{A^A)''g\ (Eq,22) 

10 where the superscript T refers to the matrix transpose, is the variance (i.e. sum of 
squares of the differences between the data and the computed fit divided by n-4). The 
diagonals of the matrix E represent the errors in each unknown so that the square root of 

the m!^ diagonal element (i;e. .y^^)is the standard error in the m'^ unknown. 

The least squares method allows for the computation of errors contained in the 
1 5 matrix A, which include both the variance (error) caused by poor data quality or random 
noise as well as the expected error resulting fi-om the data distribution. Note that the 
elements of the matrix A are functions of the offsets and azimuths in the data, which are 
then combined and used to compute the errors. In addition, since the variance caused by 
poor data quality or random noise can be computed independently, the variance caused by 
20 poor data quality or random noise, 

E^iA^A)-', (Eq.23) 

and the error resulting from the data distribution, 

E= (7 , (Eq. 24) 
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can be separated in step 104 and either one or both compared in step 106 to the AW or 
AVOA result to confirm whether or not the result has an acquisition 'footprint' - a pattern 
that is caused by the acquisition geometry. Typically, this comparison is done visually, 
although those of ordinary skill in the art would know to make the comparison 
5 mathematically. An indicia of the accuracy of the obtained results is the absence of the 

acquisition geometry (such as fold, maximum offset and minimum offset) being present in 
the obtained results. 

Errors associated with the calculation of the amplitude variation with azimuth may 
be calculated in an analogous manner to the calculation of the errors associated with the 
10 calculation of the time shift variation with azimuth, with the process being applied to Eq. 
12 rather than Eq. 19. Because the velocity variation is calculated from the travel time 
variation, the errors calculated for the time shift variations with azimuth are applicable to 
the velocity variations with azimuth. 

COMPUTING A NEW SURFACE CONSISTENT STATICS SOLUTION. 

1 5 An azimuthal velocity variation (AW) will cause offset dependent and time 

dependant statics in CMP sorted gathers. Surface consistent statics solutions for a CMP 
gather typically use a pilot trace comprising a stack which includes all offsets and azimuths 
in the gather. If AW exists, the far offsets will not align v^th the near offsets and will not 
stack coherently, thus the pilot trace will be representative of the near offset traces, which 

20 are relatively unaffected by the AW. To obtain a single static per trace, a cross- 
correlation between each trace and the pilot trace is performed over a large time window 
(for example, two seconds), typically around a target horizon. If the time delay due to 
AW changes with time, then the pilot trace will not be representative of the far offset 
traces, since the far offset traces will be stretched and squeezed relative to the near offsets 

25 and thus also relative to the pilot trace. For instance, a particular azimuth may be in the 
slow direction (resulting in a shift to later times) at one traveltime, but change to the fast 
direction (resulting in a shift to earlier times) at a later traveltime. Thus when this trace is 



-20- 



AXG-001 



correlated with the pilot trace, a poor cross-correlation may be obtained, possibly resulting 
in an incorrectly picked static for that trace. Even if a static is picked that represents the 
^average' time delay due to AW, it will appear as noise in the surface consistent statics 
computation (the "SCSC") because this computation assumes the statics are surface 
consistent, when in fact the time delays due to AW are not. 

In accordance with an embodiment of the present invention, as outlined in FIG. 7, 
an iterative process is utilized in which, in step 1 1 0, time shifts are computed as described 
above with reference to FIG. 5, which are then applied to the seismic data traces in a 
manner equivalent to a normal moveout correction, followed by step 1 12, in which surface 
consistent statics calculations known to the prior art are performed. The time shift step of 
1 10 is then repeated. This process of computing and applying the time shifts resulting 
from azimuthal velocity variations and computing the surface consistent statics is repeated 
until the process converges. This determination of whether the process has converged is 
made in step 1 14. One criterion that may be applied to determine whether the process is 
converged is whether the time shifts computed from the surface consistent statics 
computation are generally less than two milliseconds. 

At such time as the time shifts resulting from azimuthal anisotropy have been 
applied to the seismic data traces in the gather, AVO analysis as well as AVOA analysis, 
such as discussed herein with reference to FIG. 4, may be performed on the adjusted 
traces. 

The process of the invention disclosed herein is most conveniently carried out by 
writing a computer program to carry out the steps described herein on a work station or 
other conventional digital computer system of a type normally used in the industry. The 
generation of such a program may be performed by those of ordinary skill in the art based 
on the processes descried herein, FIG. 9 shows such a conventional computer system 
comprising a central processing unit 122, a display 124, an input device 126, and a plotter 
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128. The computer program for carrying out the invention will normally reside on a 
storage media (not shown) associated with the central processing imit. Such computer 
program may be transported on a CD-ROM or other storage media shown symbolically as 
storage medium 130 

The results of the calculations according this invention may be displayed with 
commercially available visualization software. Such software is well known to those of 
ordinary skill in he art and will not be further described herein. It should be appreciated 
that the resuhs of the methods of the invention can be displayed, plotted or both 

While the invention has been described and illustrated herein by reference to certain 
preferred embodiments in relation to the drawings attached hereto, various changes and 
further modifications, apart from those shown or suggested herein, may be made herein by 
those skilled in the art, vsdthout departing from the spirit of the invention, the scope of 
which is defined by the following claims. 
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